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I.  INTRODUCTION 


A  previous  report  in  this  series  [4]  described  the  simultaneous  or  subspacc  iteration 
method  for  the  partial  solution  of  the  generalized  symmetric  eigenvalue  problem  and  pro¬ 
vided  a  USA  Standard  Fortran  X3.9  -  1966  subroutine  subprogram  SIMITZ  implementing 
the  algorithm.  A  synopsis  of  that  report  was  later  published  in  a  technical  journal  [5] 
and  served  as  a  model  for  a  Fortran  X3.9  -  1978  (Fortran  77)  version  issued  in  1983  as  a 
component  of  a  proprietary  scientific  software  library  distributed  worldwide  by  Numerical 
Algorithms  Group  (NAG).  Limited,  of  Oxford.  United  Kingdom.  This  code  is  referred  to 
in  the  NAG  Fortran  Library  as  F02FJF  [6j. 

The  SIMITZ  code  as  originally  published  proved  in  early  experiments  to  be  deviant 
in  the  Fortran  77  standard.  Accordingly,  the  1966  code  has  been  rewritten  in  Fortran  77 
utilizing  all  of  the  applicable  features  of  this  standard.  The  revised  code  is  included  with 
this  report. 

Recent  experiments  with  the  Fortran  77  version  of  SIMITZ  have  been  carried  out  using 
the  CFT77  compiler  on  a  CRAY  X-MP  supercomputer.  This  compiler  offers  extensions  to 
Fortran  77  some  of  which  are  part  of  the  proposed  Fortran  8x  standard  [3].  Accordingly,  the 
Fortran  77  version  has  been  rewritten  to  incorporate  several  Fortran  8x  standard  features, 
and  that  code  is  also  included  with  this  report.  This  new  Fortran  8x  version  of  SIMITZ 
exposes  the  algorithm  for  the  first  time  to  the  full  potential  of  multiple  vector  processors 
incorporating  multiple  functional  units  featured  with  the  supercomputer  offerings  of  Cray 
Research,  Inc. 

Section  II  of  this  report  furnishes  a  brief  mathematical  description  of  subspace  iteration. 
It  also  includes  relevant  implementation  details  and  a  description  of  test  programs  to  enable 
rapid  checks  for  correct  installation.  The  author’s  recommendations  comprise  Section  III. 
Finally,  appendices  to  this  report  provide  listings  of  the  machine-readable  documentation 
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furnished  with  each  version  ol  SIMITZ.  A  5^"  diskette  which  includes  an  ASCII  file  of 
each  version  is  provided  to  qualified  requestors  of  the  present  report. 
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II  DESCRIPTION 


The  present  programs  are  implementations  of  the  subspace  iteration  algorithm  [7]  for 
calculating  the  eigenvalues  largest  in  magnitude  and  corresponding  eigenvectors  of  a  real 
matrix  symmetric  relative  to  a  prescribed  inner  product.  Let  ip(n.  u\z)  denote  an  inner 
product  in  the  space  of  real  column  n-tuples  and  let  the  real  n-square  matrix  C  satisfy 
ip(n.C u\  z)  —  ip{v .  if .  C  z) .  Then  C  is  symmetric  relative  to  ip,  and  if  the  n-square  positive 
definite  matrix  B  satisfies  ip(n.  u'.z)  —  u,TBz  then  C  is  B-syrnmetric .  The  equation  BC  — 
Cr B  characterizes  the  B-symmetry  of  C .  Given  an  optional  set  of  p  initial  approximate 
eigenvectors  of  a  real  n-square  B-symmetric  matrix  C  corresponding  to  p  eigenvalues  of 
C  largest  in  magnitude,  the  program  calculates  em  eigenvalues  and  em  corresponding 
eigenvectors  0  <  em  <  p  <  n.  to  a  precision  dependent  on  the  structure  of  C  and  on  a 
prescribed  tolerance  eps.  The  matrix  B  is  presented  to  the  program  as  an  independently 
prepared  real  function  subprogram  which  calculates  ip(n,w,z)  —  u;T Bz  given  column  n- 
vectors  u:  and  z.  The  matrix  C  is  presented  as  an  independently  prepared  subroutine 
subprogram  op(n.z.iv)  which  when  given  an  n-vector  z  computes  its  image  w  =  Cz.  The 
program  is  an  outgrowth  of  a  literal  Fortran  translation  of  the  ALGOL  procedure  ritzit  [8j 
to  which  it  is  substantially  equivalent  when  C  —  C1  and  ip(n.w.z)  =  wTz.  the  standard 
inner  product.  But  depending  on  the  choice  of  B  and  C.  the  present  program  enables  the 
direct  treatment  of  a  wide  variety  of  symmetric  eigenproblems. 

Let  A  =  At  and  B  =  BT  denote  n-square  real  matrices  and  let  o  be  real.  If  B  is 
positive  definte  then  the  matrix  C  =  B-1(A  -  oB)  is  B-symmetric,  and  the  program 
computes  eigenvalues  farthest  from  o  of  the  eigenproblem  Au  =  \Bu  and  corresponding 
eigenvectors.  Implementation  of  op(n,z,v!)  here  consists  in  providing  for  the  appropriate 
solution  for  w  of  the  linear  system  Bw  —  (A  —  oB)z.  Alternatively,  selection  of  op  to  solve 
the  system  (A-oB)u'  =  Bz  for  w  enables  the  calculation  by  simultaneous  inverse  iteration 
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of  tin  eigenvalues  nearst  to  a  and  their  eigenvectors.  Implic at  ion-  tor  large  sparse  systems 
for  which  the  Cholesky  factorization  of  B  is  impractical  are  clear. 

The  user  may  wish  to  supplement  the  following  outline  of  the  operation  of  the  program 
by  consulting  the  description  of  the  ALGOL  procedure  ritzit  in  [8]  as  well  as  a  review  of 
the  mathematical  foundations  of  simultaneous  iteration  in  [7  . 

Let  the  eigenvalues  d, . dp.dp~\ . d„  of  C  be  arranged  in  order  of  descending  ab¬ 

solute  value  and  let  Ep  denote  the  direct  sum  of  the  distinct  eigenspaces  corresponding 

to  d\ . dp.  Let  A'o  denote  an  n-by-p  matrix  having  a  /^-dimensional  column  space  not 

orthogonal  relative  to  ip  to  any  eigenvector  in  Ep.  Simultaneous  iteration  is  based  on  the 
observation  that  if  idpj  >  idp_i|.  the  columns  of  the  matrix  Ah~m  =  CmAh  tend  to  a  basis 
of  Ep  as  ks  =  k  +  m  increases.  But  in  practice  all  of  the  columns  of  A' *  tend  toward  the 
eigenspace  Ex  causing  loss  of  information  concerning  the  residual  eigenvectors.  To  counter 
this  tendency,  set 

A'*-m  =  CmXkRk}m  (1) 

where  the  p-square  upper  triangular  matrix  Rk~m  >s  constructed  together  with  X^^m  by 
the  Gram-Schmidt  process  to  render  the  columns  of  A\-  orthonormal  relative  to  ip.  Now 
the  ith  column  vector  of  A”/-  converges  to  the  ith  eigenvector  of  C  at  a  rate  proportional 
to  max2<,<p(|d|/d,_i  j.  \dt .  x/d,\).  Clearly  this  convergence  will  be  delayed  in  the  presence 
of  eigenvalue  clustering.  But  if  \dp‘  —  |dp_i|  is  not  too  small,  the  column  space  of  Ah  will 
contain  a  good  approximation  to  the  ith  eigenvector  even  when  k  is  small. 

In  order  to  recover  this  approximation,  a  modified  Rayleigh-Ritz  process  is  employed. 
Let  Qk  denote  an  orthogonal  matrix  which  diagonalizes  the  p-square  symmetric  matrix 
R^Rf .  Then  the  itr>  column  vector  of 

.A,.,  -  r\kR,\Qk  ,  (2) 
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coverges  to  the  i!h  eigenvector  of  C  at  a  rate  proportional  to  dp .  ,  jd,  while  the  entries  of 
the  diagonal  matrix  computer  with  Qk  and  properly  ordered  offer  close  approximations  to 
d\.  .  ..dp.  The  true  signed  eigenvalues  need  only  oe  computed  at  termination  by  diagonal¬ 
izing  the  leading  [p  —  1  l-square  principal  suhmatrix  of  X[.  B<  -\\.  the  eigenproblem  for  (' 
projected  on  Efl  ,  relative  to  ip. 

The  program  determines  a  strategy  for  employing  the  devices  (1)  and  (2)  based  on  the 
distribution  of  the  leading  p  eigenvalues  of  ('  upon  which  the  convergence  rate  ultimately 
depends.  The  selection  of  values  w  in  (1)  is  particularly  important  in  this  regard  in  that 
CmXi(  is  replaced  by  the  rn,h  Chebychov  polynomial  on  the  inter.nl  ;  —  e.e’  evaluated  by  a 
special  3-term  recurrence  relation  and  permitting  accelerated  convergence  when  values  of 
rn  are  continually  large:  f  is  the  current  value  of  dp.  As  a  result  the  convergence  quotient 
lies  between  \dp/dern\  and  exp( -arccosh  \demjdp\).  It  is  nearer  to  the  first  value  if  \d\/dtrn\ 
is  large  and  nearer  to  the  second  if  the  latter  quotient  is  close  to  one. 

As  the  iteration  proceeds  through  a  maximum  of  j/rml  iteration  steps  --  km  is  a  pro¬ 
gram  parameter  —  acceptance  tests  for  the  eigenvalues  and  eigenvectors  are  conducted 
following  each  of  the  Rayleigh-Ritz  steps  (2).  As  soon  as  the  relative  increase  of  |o'/)-1( 
is  smaller  than  tps/10.  then  d) is  accepted  and  h .  the  number  of  previously  accepted 
eigenvalues,  is  increased  by  one.  Eigenvectors  are  accepted  in  groups  of  one  or  more  corre¬ 
sponding  to  clusters  of  accepted  eigenvalues  nearly  equal  in  magnitude.  If  g  eigenvectors 
have  already  been  accepted,  let  dg^}....d(  denote  such  a  cluster.  For  all  j,  g  +  1  <  j <  t, 
denote  by  y}  the  projection  relative  to  ip  of  the  image  Ci;  of  the  jth  column  x7  of  X^  on 

the  linear  closure  of  Xj . x/.  Set  /,  =  max7  | \Cx}  —  y}\\/\\C x.j\\  for  i  -  g  + 1, ...,  f  where  the 

indicated  norm  is  the  Euclidean  norm  or  2-norm  relative  to  ip.  If  \d(\f ij {\d(\  —  e)  is  smaller 
than  tps  then  all  the  x} .  j  =  g  - 1 1 .....  f,  are  accepted  as  eigenvectors  and  g  is  increased  to  L 
The  error  quantities  /,  are  systematically  discounted  in  accordance  with  the  convergence 


:> r t >;•<”* i t  ~  of  t:.»  . t . c. ■  > r : ' : ; : i i  *«•  permit  nmvcr^'iici1  in  the  presence  of  excessive  round-off 
error  or  >m><  ’  parameter  t p.<  is  prescribed  unrealistically  small.  Having  determine  g 
eig,  ;•  vector',  'in  it«  ration  continues  with  p  -  y  remaining  columns  of  A”*  until  either  ttn 
eigetivecti  r»  i.iivi  been  calculated  or  A-?m  has  been  exceeded.  The  program  may  reduce 
t  tt>  if  it  defects  <  ;:iit  r  no  progress  in  convergence  of  eigenvectors  corresponding  to  smaller 
eigenvalue'  or  hick  <*f  stability  m  the  behavior  of  larger  eigenvalues. 

The  Fortran  programs  presented  here  which  implement  the  algorithm  described  above 
differ  primarily  it.  their  treatment  of  loci  storage  required  by  SIMITZ  and  in  the  parame¬ 
ters  associated  with  the  procedure  op[  n.  z.  w).  The  Fortran  77  version  requires  the  calling 
program  to  reserve  max|p'.2n)  —  3 p  consecutive  locations  of  temporary  storage  for  use  by 
SIMITZ.  The  CFT77  code  tends  to  trade  memory  for  speed.  There  temporary  storage  is 
accommodated  in  Fortran  Sx  "automatic  arrays"  which  total  2 np  +  p2  -+•  4p  locations  and 
do  not  involve  the  calling  program.  Both  versions  are  adapted  to  operate  in  the  software 
environment  furnished  by  the  SLATEC  Ma’hematical  Subprogram  Library  [9],  Version 
3.1.  Calling  sequences  arc  identical  and  agree  with  that  of  the  original  Fortran  6G  code. 

The  role  of  the  procedure  op{n.z.  w)  in  the  CFT77  code  for  the  CRAY  differs  slightly 
from  that  in  the  Fortran  66  and  Fortran  77  versions.  In  the  former  code  z  and  ir  identify 
matrices  of  column  n-vectors.  and  op  must  compute  the  image  tr  of  r  under  multiplication 
by  the  matrix  ('.  Hence,  op  requires  two  additional  parameters  which  identify  (1)  the 
extent  of  the  array  containing  the  matrix  z  and  (2)  the  number  of  column  vectors  of  the 
matrices  identified  bv  .:  and  u\  Details  may  be  found  in  Appendix  B.  A  view  of  what  can 
be  accomplished  when  op  implements  matrix  multiplication  by  a  full  matrix  C  may  be 
inferred  from  a  paper  of  Bailey  1  .  More  complicated  computations  involving  the  matrices 
('  and  c  require  a  more  elaborate  and  careful  design  of  the  procedure  op.  A  survey  by 
Duff  2.  conveys  the  flavor  of  several  of  these  design  issue'  when  the  targft  computer  is  a 

c. 


multiple  processor  CRAY. 

The  Fortrai;  77  version  of  SIMITZ  provided  with  the  present  report  was  installed  and 
revised  on  a  CDC  CYBER  180-S55  processor  using  the  NOS  2.5  operating  system  and  the 
FTN5  Fortran  compiler.  Subsequently,  installation  was  carried  f'ut  on  a  CRAY  X-MP/12 
processor  with  the  COS  1.15  operating  system  and  CFT  1.15  Fortran  compiler.  Revisions 
adhering  to  the  Fortran  8x  standard  were  completed  with  the  CRAY  CFT77  compiler. 
Version  1.3. 

The  revision  process  was  accomplished  in  part  with  the  aid  of  an  executable  Fortran 
program  TESTD  (liven  as  input  the  integer  values  n.p.  and  km  together  with  a  real  value 
e/)s.  TESTD  generates  a  pair  of  random  real  n-square  diagonal  matrices  A  and  B .  It  then 
exercise?  SIMITZ  on  the  eigenvalue  problem  AX  -  BXD  =  O  where  D  =  diag(d\ . dp). 
For  each  value  em\in)  =  ern.ern  =  1.  TESTD  calls  SIMITZ  and  records  the 

number  ern(r)  of  eigenpairs  successfully  calculated  within  the  tolerance  eps ,  the  maximum 

fmaT(r)  of  the  error  quantities  /, ,  i  —  1 . ern(r).  and  the  total  number  Es(r)  of  iterations 

expended  by  SIMITZ  in  each  call.  Each  exercise  is  performed  twice,  once  for  random 
(r)  initial  eigenvectors  A”0  and  once  when  initial  eigenvectors  are  Lanczos  (£)  vectors  [7], 
TESTD  finally  lists  the  true  eigenvalues,  which  are  the  diagonal  entries  c,  of  B  lA.  and 
their  computed  counterparts  dt.  the  corresponding  error  quantities  /,  and  the  norms  e, 

in  the  standard  Euclidean  metric  of  the  residuals  (.4  -  dtB)X.i  =  1 . em(f).  Figure  1 

shows  a  typical  output  listing  from  the  executable  program  TESTD. 
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Fig.  1.  Output  List(CFT77) 


N  =  100 

P  =  10 

KM 

=  210 

EPS  =  1 

.  25E-06 

EM(IN) 

EMIL) 

FMAX(L) 

KS  (L) 

KS  (R) 

FMAX(R) 

EM(R) 

1 

1 

2.24E-08 

17 

30 

4. 13E-09 

1 

2 

2 

1 .05E-07 

20 

30 

1.87E-09 

2 

3 

3 

1 .75E-07 

43 

50 

1 . 55E-08 

3 

4 

4 

1 .48E-07 

56 

56 

8.45E-08 

4 

5 

5 

2.49E-07 

75 

75 

5.36E-08 

5 

6 

6 

1 . 19E-07 

94 

75 

1 . 14E-08 

6 

7 

•7 

2.61E-07 

113 

138 

3.51E-08 

6 

8 

8 

2.36E-07 

113 

138 

7.71E-08 

6 

9 

6 

1 .64E-07 

138 

138 

8. 10E-08 

6 

TOTAL  KS (L) 

=  669 

730  = 

TOTAL  KS  (R) 

I 

C  ( I ) 

D  ( I ) 

F  (I) 

E(I) 

1  -1 . 545762930 1304E+00  -1 . 545762930 1305E+00  7.29E-08  9.84E-08 

2  1 .46043625 16858E+00  1 . 46043625 16858E+00  1.64E-07  2.28E-07 

3  -1. 1321574055138E+00  - 1 . 132 1574055 138E+00  5.27E-08  5.82E-08 

4  1 . 1042303464 189E+00  1 . 1042303464 190E+00  9.61E-10  9.55E-10 

1 . 0432976702273E+00  1 . 0432976702274E+00  1 . 37E-08  1 . 13E-08 

-1 .0257765390463E+00  - 1 . 0257765390463E+00  5.10E-08  3.60E-08 
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III.  recommendation  - 


The  computist  having  an  application  requiring  partial  solution  of  the  large  sparse 
symmetric  eigen t  aiue  problem  now  has  several  options.  If  he  has  access  to  the  NAG  Fortran 
Library.  Mark  1 1 .  or  to  a  later  mark,  he  would  do  well  to  consider  use  of  F02FJF  from  that 
library.  The  NAG  version  of  SIMITZ  has  benefited  by  extensive  independent  testing  by  an 
internationally  known  technical  staff  as  well  as  from  field  testing  by  a  worldwide  clientele. 
Moreover,  a  clever  use  of  the  singular  value  decomposition  avoids  the  computation  of  the 
squares  of  the  desired  eigenvalues  thus  enhancing  robustness  ot  the  code  when  the  range 
of  Fortran  real  values  may  be  limited  as  on  some  microcomputers.  This  writer  was  unable 
to  match  the  accuracy  of  the  NAG  auxiliary  routines  for  singular  value  decomposition 
of  an  upper  triangular  matrix  using  the  LINPACK  code  SSVDC  found  in  the  SLATEC 
collection.  F02FJF  may  run  a  bit  slower  on  some  processors  than  its  Fortran  77  SIMITZ 
counterpart  owing  to  extensive  subprog-am  calls  to  NAG  auxiliary  routines.  This  speed 
differentia!  will  decrease  as  the  size  n  of  the  matrix  increases. 

If  one  desires  to  use  codes  which  are  all  in  the  public  domain,  he  may  employ  either 
the  Fortran  77  SIMITZ  or  the  Fortran  8x  version  for  CRAY  processors.  Accuracy  and 
robustness  are  almost  identical  owing  to  adherence  faithfully  to  the  ANSI  standard  and  to 
the  use  of  the  SLATEC  auxiliary  programs. 

Whichever  version  one  uses,  he  should  need  the  following  advice:  Gain  experience  with 
SIMITZ  on  problems  of  moderate  size  before  investing  extensive  computer  resources  on 
large  problems  for  which  the  codes  are  ultimately  intended. 


9 


REFERENCES 


1.  David  H.  Bailey.  Extra  high  speed  matrix  multiplication  on  the  CRAY-2.  SIAM  J.  Sci. 
Stat.  Comput.  9  (1988).  603-607. 

2.  Iain  S.  Duff.  The  influence  of  vector  and  parallel  processors  on  numerical  analysis, 
in  The  State  of  the  Art  in  Numerical  Analysis,  A.  Iserles  and  M.J.D.  Powell,  eds., 
Clarendon.  Oxford,  1987. 

3.  Michael  Metcalf  and  John  Reid.  Fortran  8x  Explained.  Clarendon,  Oxford.  1987. 

4.  Paul  J.  Nikolai.  Solving  the  real  generalized  symmetric  eigenproblem  by  simultane¬ 
ous  iteration.  Technical  Report  AFFDL-TR-76-118  (1976),  distributed  by  National 
Technical  Information  Service,  U.S.  Department  of  Commerce.  5285  Port  Royal  Road, 
Springfield.  VA  22161. 

5.  Paul  J.  Nikolai,  ALGORITHM  538,  Eigenvectors  and  eigenvalues  of  real  generalized 
symmetric  matrices  by  simultaneous  iteration,  ACM  Trans.  Math.  Software  5  (1979), 
118-125. 

6.  Numerical  Algorithms  Group,  Ltd.,  The  NAG  Fortran  Library  Manual.  Mark  11,  Nu¬ 
merical  Algorithms  Group.  Downers  Grove,  1983. 

7.  Beresford  N.  Parlett,  The  Symmetric  Eigenvalue  Problem,  Prentice-Hall,  Englewood 
Cliffs,  1980. 

8.  H.  Rutishauser,  Simultaneous  iteration  method  for  symmetric  matrices,  Numer.  Math. 
16  (1970).  205-223. 

9.  Walter  H.  Vandevender,  The  SLATEC  Mathematical  Subprogram  Library.  Version 
2.0.  Sandia  Report  SAND-84-0281,  Sandia  National  Laboratories.  Albuquerque,  1984, 
distributed  by  National  Technical  Information  Service,  U.S.  Department  of  Commerce. 
5285  Port  Royal  Road.  Springfield,  VA  22161. 


10 


APPENDIX  A 


Appendix  A  provides  a  listing  in  SLATEC  format  of  the  Fortran  77  documentation  of 

SIMITZ. 
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C***BEGIN  PROLOGUE  SIMITZ 
C*#*DATE  BITTEN  750815  (YYNMDD) 

C**#REVISION  DATE  881022  (YYMVDD) 

C  *  *  *  C  AT AGORY  NO.  F2C2 ,  F2C9 ,  F2D 

C  *  #  *  KEYWORDS  EIGENVALUES, EIGENVECTORS, SUBSPACE  ITERATION 
C##*AUTHOR  NIKOLAI,  PAUL  J. 

C  U.S.  AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 

C  WRIGHT-PATTERSON  AFB ,  OH  45433 

C#»#PURPOSE  GIVEN  AS  OPTIONAL  INPUT  A  SET  OF  P  INITIAL  APPROXIMATE 
C  EIGENVECTORS  OF  A  REAL  N- SQUARE  SYMMETRIC  MATRIX  A  CORRES- 

C  PONDING  TO  P  EIGENVALUES  LARGEST  IN  MAGNITUDE,  SIMITZ  COM- 

C  PUTES  EM  EIGENVALUES  LARGEST  IN  MAGNITUDE  AND  EM  CORRES- 

C  PONDING  EIGENVECTORS  TO  A  PRECISION  DEPENDING  ON  THE  STRUC- 

C  TORE  OF  A  AND  ON  A  PRESCRIBED  TOLERANCE  EPS.  THIS  VERSION 

C  OF  SIMITZ  IS  AN  ANSI  X3.9-1978  FORTRAN  VERSION  OF  (1) . 

C**#DESCRIPTION 
C 

C  CONTROL 
C 

C  DIMENSION  X(LDX,P) ,  D(P) ,  WK(K) 

C  INTEGER  P,  EM 

C  REAL  IP 

C  EXTERNAL  IP,  INF,  OP 

C 
C 
C 

C  CALL  SIMITZ (N,  P,  KM,  EPS,  IP,  OP,  INF,  EM,  X,  LDX,  D,  WK) 

C 

C  WHERE 

C  N  IS  AN  INTEGER  INPUT  VARIABLE,  THE  ORDER  OF  THE  MATRIX  A. 

CP  IS  AN  INTEGER  INPUT  VARIABLE,  THE  NUMBER  OF  SIMULTANEOUS 

C  ITERATION  VECTORS . 

C  KM  AS  AN  INTEGER  INPUT  VARIABLE  IS  IN  MAGNITUDE  THE  MAXIMUM 

C  NUMBER  OF  ITERATION  STEPS  TO  BE  EXECUTED.  IF  KM  IDENTIFIES 

C  A  NEGATIVE  VALUE  THEN  P  INITIAL  APPROXIMATE  EIGENVECTORS 

C  ARE  ASSUMED  TO  BE  PRESENT  IN  THE  ARRAY  X.  OTHERWISE  SIMITZ 

C  SUPPLIES  RANDOM  INITIAL  EIGENVECTORS. 

C  KM  AS  AN  INTEGER  OUTPUT  VARIABLE  IDENTIFIES  THE  NUMBER  KS  OF 
C  ITERATION  STEPS  FINALLY  USED  IN  THE  CALCULATION  OF  EM 

C  EIGENVECTORS . 

C  EPS  IS  A  REAL  INPUT  VARIABLE,  THE  TOLERANCE  FOR  ACCEPTING 
C  EIGENVECTORS.  AS  SOON  AS  SUCCESSIVE  ITERATES  OF  THE  RITZ 

C  VALUES  ABS(D(H+1) )  DIFFER  BY  LESS  THAN  ABS(D(H+1) ) *EPS/10.0 

C  THEN  D(H+1)  IS  ACCEPTED  AS  AN  EIGENVALUE  AND  H,  THE  NUMBER 

C  OF  PREVIOUSLY  ACCEPTED  EIGENVALUES,  IS  INCREASED  BY  1 .  AS 

C  SOON  AS  THE  ERROR  QUANTITIES  F(I) ,  NORMS  OF  THE  RESIDUALS. 

C  SATISFY  D (I) *F(I) / (D (I )  -  D(P)1  .LT.  EPS,  THEN  G,  THE  NUM- 

C  BER  OF  ALREADY  ACCEPTED  RITZ  VECTORS,  IS  INCREASED  TO 

C  G  +  L,  I  =  G  +  1 . L.  THE  F(I)  ARE  DISCOUNTED  WITH 

C  SUCCESSIVE  ITERATIONS  TO  FORCE  CONVERGENCE  IN  CASE  OF  UN- 

C  FORTUNATE  CHOICE  OF  PARAMETERS.  IF  M  SIGNIFICANT  DIGITS 

C  OF  ACCURACY  ARE  REQUIRED  OF  THE  EIGENVALUES,  THEN  SET 

C  EPS  EQUAL  TO  10.0**(-M)  AS  A  GENERAL  RULE. 

C  IP  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A  FORTRAN  COM- 
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C  PATIBLE  REAL,  FUNCTION  SUBPROGRAM  OF  THE  FORM  IP(N,  2,  W) 

C  WHICH  MUST  RETURN  THE  INNER  PRODUCT  OF  THE  VECTORS  IDENTI- 

C  FT  ED  BY  THE  N- ARRAYS  Z  AND  W.  WHEN  A  IS  SYMMETRIC,  IP 

C  SHOULD  RETURN  THE  STANDARD  INNER  PRODUCT  ( Z- TRANSPOSED )  W. 

C  WHEN  A  IS  SYMMETRIC  RELATIVE  TO  A  GENERAL  INNER  PRODUCT 

0  ( Z  ■  TRANS ;  'OSEI  > '  BW .  B  POSITIVE  DEFINITE.  THEN  IP  MUST  RETURN 

0  THIS  INNER  PRODUCT . 

C  OP  IS  AN  EXTERNAL  INPUT  VARIABLE.  THE  NAME  OF  A  FORTRAN  COM- 
C  PATIBLE  SUBROUTINE  SUBPROGRAM  OF  THE  FORM  OP  IN,  Z.  W) 

C  WHICH  MUST  CALCULATE  THE  IMAGE  W  OF  THE  VECTOR  IDENTIFIED 

C  BY  THE  N- ARRAY  Z  UNDER  THE  N-SQUARE  MATRIX  A  WITHOUT  OVER- 

C  WRITING  2. 

C  INF  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A  FORTRAN  COM¬ 
'S  PATIBLE  SUBROUTINE  SUBPROGRAM  WHICH  MAY  BE  USED  FOR 

C  OBTAINING  INFORMATION  OR  TO  EXERT  CONTROL  DURING  EXECUTION 

C  OF  SIMIT2.  INF  HAS  THE  FORM  INFCKS,  G,  H,  F)  WHERE 

C  KS  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  OF  THE  NEXT 

C  ITERATION  STEP. 

C  G  IS  AN  INTEGER  OUTPUT  VARIABLE.  THE  NUMBER  OF  ALREADY 

C  ACCEPTED  EIGENVECTORS. 

C  H  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  OF  ALREADY 

C  ACCEPTED  EIGENVALUES. 

C  F  IS  A  REAL  OUTPUT  VARIABLE  P-ARRAY,  ERROR  QUANTITIES 

C  MEASURING  RESPECTIVELY  THE  STATE  OF  CONVERGENCE  OF 

C  THE  P  SIMULTANEOUS  ITERATION  VECTORS.  EACH  ELEMENT  OF 

C  THE  ARRAY  F  IS  INITIALLY  SET  BY  SIMITZ  TO  THE  VALUE  4.0. 

C  EM  AS  AN  INTEGER  INPUT  VARIABLE  IS  THE  NUMBER  OF  EIGENVALUES 
C  TO  BE  COMPUTED ,  0  .LT.  EM  .LT.  P  .LE.  N  .LE.  LDX. 

C  EM  AS  AN  INTEGER  OUTPUT  VARIABLE  IS  THE  NUMBER  OF  EIGENVECTORS 
C  COMPUTED  THROUGH  KM  ITERATION  STEPS. 

C  X  AS  A  REAL  N-BY-P  INPUT  ARRAY  IS  A  SET  OF  P  OPTIONAL  INITIAL 

C  APPROXIMATE  EIGENVECTORS  X(I,1),  _ X(I,P>,  I  =  1 . 

C  N,  INTERPRETED  BY  SIMITZ  IF  KM  IS  NEGATIVE. 

C  X  AS  A  REAL  N-BY-P  OUTPUT  ARRAY  IS  A  SET  OF  EM  EIGENVECTORS 

C  X  ( I  ,  1 ) . X  ( I  ,  EM)  ,  I  =  1  ,  ....  N ,  COMPUTED  THROUGH 

C  ABS (KM)  ITERATION  STEPS  WITH  THE  REMAINDER  OF  X  CONSISTING 

C  OF  P  -  EM  APPROXIMATE  EIGENVECTORS.  THE  P-SQUARE  MATRIX 

C  WOSE  (J,  L)  ENTRY  IS  IP(N,  X(l.J),  X(1,L))  IS  THE  IDEN- 

C  TITY  MATRIX;  THAT  IS,  THE  EIGENVECTORS  OF  A  ARE  ORTHO- 

C  NORMAL  RELATIVE  TO  THE  INNER  PRODUCT  IP. 

C  LDX  IS  AN  INTEGER  INPUT  VARIABLE  V«ICH  IDENTIFIES  THE  LEADING 
C  DIMENSION  IN  THE  CALLING  PROGRAM  OF  THE  ARRAY  X. 

CD  IS  A  REAL  OUTPUT  P-ARRAY  OF  WHICH  D(l) ,  ....  D (EM)  ARE  THE 

C  EIGENVALUES  OF  A  LARGEST  IN  MAGNITUDE  IN  DECREASING  ORDER 

C  CORRESPONDING  TO  THE  COMPUTED  EIGENVECTORS  X(I , 1) ,  .... 

C  X ( I , EM)  ,1  =  1 . N.  D(EM+1) . D(P-l)  CONTAIN 

C  APPROXIMATIONS  TO  PROGRESSIVELY  SMALLER  SUCH  EIGENVALUES. 

C  D(P)  CONTAINS  THE  MDST  RECENTLY  COMPUTED  VALUE  OF  E,  WHERE 

C  THE  INTERVAL  (-E,  E)  IS  THE  INTERVAL  OVER  WHICH  THE 

C  CHEBYSHEV  ACCELERATION  WAS  PERFORMED. 

C  WK  THE  INITIAL  LOCATION  OF  AT  LEAST  MAX(P**2,  2*N)  +  3»P  =  K 
C  CONSECUTIVE  STORAGE  LOCATIONS  WHICH  MAY  NOT  BE  OVER- 

C  WRITTEN  WHILE  SIMITZ  IS  IN  EXECUTION. 

C 


C  OTHER  PROGRAMMING  INFORMATION 
C 

C  THIS  VERSION  OF  SIMITZ  IS  DESIGNED  TO  OPERATE  IN  THE  SOFTWARE 
C  ENVIRONMENT  FURNISHED  BY  THE  SLATEC  COLLECTION,  VERSION  3.1. 

C 

C  THE  PERFORMANCE  OF  SIMITZ  IS  STRONGLY  DEPENDENT  UPON  THE  CHOICE 
C  OF  INPUT  PARAMETERS  AND  UPON  THE  CAREFUL  PREPARATION  OF  THE 
C  SUBPROGRAMS  IV  AND  OP.  THE  USER  SHOULD  CONSIDER  USING  HIS  OWN 
C  ACTIVE  SUBROUTINE  INF  TO  MONITOR  PROGRESS  OF  SIMITZ  RELATIVE  TO 
C  HIS  CHOICE  OF  INPUT  PARAMETERS  IF  NO  INFORMATION  IS  OTHERWISE 
C  AVAILABLE  CONCERNING  THE  LOCATIONS  OF  THE  RELEVANT  EIGENVALUES. 

C  RECALL  THAT  SIMITZ  MAY  BE  REENTERED  WITH  KM  .LT.  0  WITHOUT  LOSS 
C  OF  INFORMATION  TO  PERMIT  CONSERVATIVE  INITIAL  CHOICES  OF 
C  ABS  (KM)  ,  EPS  AND  P. 

C 

C  OTHER  PROGRAMS  REQUIRED 
C 

C  FUNCTION  RAND 

C  RETURNS  UNIFORMLY  DISTRIBUTED  RANDOM  NUMBERS  ON  THE  OPEN 

C  INTERVAL  (0,  1). 

C  SUBROUTINE  TRED2 

C  IS  THE  El SPA CK  (4)  PROGRAM  WHICH  COMPUTES  A  HOUSEHOLDER 

C  TRIDIAGONAL  FORM  OF  A  REAL  SYMMETRIC  MATRIX. 

C  SUBROUTINE  IMTQL2 

C  IS  THE  EISPACK  PROGRAM  WHICH  COMPUTES  THE  EIGENVALUES  AND 

C  ORTHONORMAL  EIGENVECTORS  OF  A  SYMMETRIC  TRIDIAGONAL  MATRIX. 

C  SUBROUTINE  XERRWV 

C  PROCESSES  AN  ERROR  (DIAGNOSTIC)  MESSAGE. 

C  FUNCTION  R1MACH 

C  RETURNS  SINGLE  PRECISION  MACHINE  DEPENDENT  CONSTANTS. 

C  FUNCTION  IP 

C  IS  DESCRIBED  ABOVE. 

C  SUBROUTINE  OP 

C  IS  DESCRIBED  ABOVE. 

C  SUBROUTINE  INF 

C  IS  DESCRIBED  ABOVE. 

C 

C  METHOD 
C 

C  SIMITZ  REPRESENTS  RESULTS  OF  EXTENSIVE  MODIFICATIONS  AND  TESTS 

C  OF  SUBROUTINE  SIMITZ  (1) ,  A  FORTRAN  66  TRANSLATION  OF  THE 

C  ALGOL  60  PROCEDURE  RITZIT  (3) .  THE  BASIC  RUT I SHAUSER- RE I NSCH 

C  ALGORITHM  IS  PRESERVED. 

C*  PREFERENCES  (1)  PAUL  J.  NIKOLAI,  ALGORITHM  538  -  EIGENVECTORS  AND 
C  EIGENVALUES  OF  REAL  GENERALIZED  SYMMETRIC  MATRICES  BY 

C  SIMULTANEOUS  ITERATION,  ACM  TRANS.  MATH.  SOFTWARE  5 

C  (1979) ,  118-125. 

C  (2)  BERESFORD  N.  PARLETT ,  THE  SYMMETRIC  EIGENVALUE 

C  PROBLEM,  PRENT ICE- HALL ,  ENGLEWOOD  CLIFFS,  1930. 

C  (3)  HEINZ  RUT I SHAUSER,  SIMULTANEOUS  ITERATION  METHOD  FOR 

C  SYMMETRIC  MATRICES,  NUMER.  MATH.  16(1970),  205-223. 

C  (4)  B.T.  SMITH  ET  AL ,  MATRIX  EIGENSYSTEM  ROUTINES  - 

C  EISPACK  GUIDE,  2-ND  ED. .LECTURE  NOTES  IN  COMPUTER 

C  SCIENCE  6.  SPRINGER-VERLAG.  NEW  YORK.  1976. 
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C***ROUTINES  CALLED  IMTQL2 ,R1MACH,RAND,THED2 .XERRWV 
C*»*END  PROLOGUE 


15 


APPENDIX  B 


Appendix  B  provides  a  listing  in  SLATEC  format  of  the  Fortran  8x  documentation  of 
SIMITZ. 


1G 


C#»#BEGIN  PROLOGUE  SIMITZ 
C***DATE  WRITTEN  750815  tYYMZDD) 

C***REVISION  DATE  881115  (YYNMDD) 

C***CATAGORY  NO.  F2C2 ,  F2C9 ,  F2D 

C*** KEYWORDS  E I GENVALUES , E I GENVECTORS , SUBSPACE  I TERAT I ON 
C**#  AUTHOR  NIKOLAI,  PAUL  J . 

C  US.  AIR  FORCE  WRIGHT  AERONAUTICAL  LABORATORIES 

C  WRIGHT -PATTERSON  AFB ,  OH  45433 

C** » PURPOSE  GIVEN  AS  OPTIONAL  INPUT  A  SET  OF  P  INITIAL  APPROXIMATE 
C  EIGENVECTORS  OF  A  REAL  N-SQUARE  SYMMETRIC  MATRIX  A  CORRES- 

C  PONDING  TO  P  EIGENVALUES  LARGEST  IN  MAGNITUDE,  SIMITZ  COM- 

C  PUTES  EM  EIGENVALUES  LARGEST  IN  MAGNITUDE  AND  EM  CORRES- 

C  PONDING  EIGENVECTORS  TO  A  PRECISION  DEPENDING  ON  THE  STRUC- 

C  TURE  OF  A  AND  ON  A  PRESCRIBED  TOLERANCE  EPS.  THIS  VERSION 

C  OF  SIMITZ  IS  A  CRAY  CFT77  FORTRAN  VERSION  OF  (1) . 

C*  ^DESCRIPTION 
C 

C  CONTROL 
C 

C  REAL  XfLDX, P) ,  DIP) 

C  INTEGER  P,  EM 

C  REAL  IP 

C  EXTERNAL  IP,  INF,  OP 

C 
C 
C 

C  CALL  SIMITZ (N,  P,  KM,  EPS,  IP,  OP,  INF,  EM,  X,  LDX,  D,  WK) 

C 

C  WHERE 

C  N  IS  A N  INTEGER  INPUT  VARIABLE.  THE  ORDER  OF  THE  MATRIX  A. 

CP  IS  AN  INTEGER  INPUT  VARIABLE,  THE  NUMBER  OF  SIMULTANEOUS 

C  ITERATION  VECTORS. 

C  KM  AS  AN  INTEGER  INPUT  VARIABLE  IS  IN  MAGNITUDE  THE  MAXIMUM 

C  NUMBER  OF  ITERATION  STEPS  TO  BE  EXECUTED.  IF  KM  IDENTIFIES 

C  A  NEGATIVE  VALUE  THEN  P  INITIAL  APPROXIMATE  EIGENVECTORS 

C  ARE  ASSUMED  TO  BE  PRESENT  IN  THE  ARRAY  X.  OTHERWISE  SIMITZ 

C  SUPPLIES  RANDOM  INITIAL  EIGENVECTORS. 

C  KM  AS  AN  INTEGER  OUTPUT  VARIABLE  IDENTIFIES  THE  NUMBER  K S  OF 
C  ITERATION  STEPS  FINALLY  USED  IN  THE  CALCULATION  OF  EM 

C  EIGENVECTORS . 

C  EPS  IS  A  REAL  INPUT  VARIABLE,  THE  TOLERANCE  FOR  ACCEPTING 
C  EIGENVECTORS.  AS  SOON  AS  SUCCESSIVE  ITERATES  OF  THE  RITZ 

C  VALUES  ABS(D(H+1) )  DIFFER  BY  LESS  THAN  ABS (D (H+ 1 ) ) #EPS/ 10 . 0 

C  THEN  DCK+l)  IS  ACCEPTED  AS  AN  EIGENVALUE  AND  H,  THE  NUMBER 

C  OF  PREVIOUSLY  ACCEPTED  EIGENVALUES.  IS  INCREASED  BY  1 .  AS 

C  SOON  AS  THE  ERROR  QUANTITIES  F(I) ,  NORMS  OF  THE  RESIDUALS, 

C  SATISFY  D(I)*F(I)/(D(I)  -  DIP))  . LT.  EPS,  THEN  G,  THE  NUM- 

C  BER  OF  ALREADY  ACCEPTED  RITZ  VECTORS,  IS  INCREASED  TO 

C  G  +  L,  I  =  G  +  1,  _  L.  THE  F(I)  ARE  DISCOUNTED  WITH 

C  SUCCESSIVE  ITERATIONS  TO  FORCE  CONVERGENCE  IN  CASE  OF  UN- 

C  FORTUNATE  CHOICE  OF  PARAMETERS.  IF  M  SIGNIFICANT  DIGITS 

C  OF  ACCURACY  ARE  REQUIRED  OF  THE  EIGENVALUES,  THEN  SET 

C  EPS  EQUAL  TO  10.0** (-M)  AS  A  GENERAL  RULE. 

C  IP  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A  FORTRAN  COM- 
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a  a  a  a  a  a 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r, 

c 

c 

c 

c 

c 

c 

n 

c 

c 


c 


? AT  I DLL  1- UNCTION  SUBPROGRAM  OF  THE  FORM  IP(N,  z,  W ) 

WICH  ML-  ."I  H.-T  LPJN  THE  INNER  PRODUCT  OK  THE  VECTORS  IDENTI¬ 
FIED  BY  7HK  N-  ARRAYS  Z  AND  W.  WEN  A  IS  SYMMETRIC ,  IP 
SHOULD  RETURN  I  HE  STANDARD  INNER  PRODUCT  ( Z -TRANSPOSED) W. 

WEN  A  OYNWETR  I C  RELATIVE  TO  A  GENERAL  INNER  PRODUCT 

iZ-7RANU\UIC'LW.  B  POSITIVE  DEFINITE.  THEN  IP  MUST  RETURN 
THIS  INNEF  rKO-DL’CT . 

OP  IF  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A  FORTRAN  COM¬ 
PATIBLE  SIR  ROUTINE  SUBPROGRAM  OF  THE  FORM  OP(N.  Z,  LDZ ,  W,  M) 
WHICH  MUST  CALCULATE  THE  IMAGE  W  OF  THE  N-BY-M  MATRIX  I DENT I - 
FIED  BY  THE  LDZ-BY-M  ARRAY  Z  UNDER  THE  N- SQUARE  MATRIX  A  WITH¬ 
OUT  OVERWRITING  Z.  SUBROUTINE  OP(N,  Z,  LDZ,  W,  M)  MUST  INCLUDE 
SPECIFICATION  STATEMENTS  EQUIVALENT  TO:  REAL  Z(LDZ,  M)  ,  W(N,  M) 

INF  IS  AN  EXTERNAL  INPUT  VARIABLE,  THE  NAME  OF  A  FORTRAN  COM¬ 
PATIBLE  SUBROUTINE  SUBPROGRAM  WICH  MAY  BE  USED  FOR 
OBTAINING  INFORMATION  OR  TO  EXERT  CONTROL  DURING  EXECUTION 
OF  S I MI TZ .  INF  HAS  THE  FORM  INF (KS,  G,  H,  F)  WHERE 

KS  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  OF  THE  NEXT 
I  '^'RAT I  ON  STEP . 

a  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUNBER  OF  ALREADY 
ACCEPTED  EIGENVECTORS. 

H  IS  AN  INTEGER  OUTPUT  VARIABLE,  THE  NUMBER  OF  ALREADY 
ACCEPTED  EIGENVALUES. 

F  IS  A  REAL  OUTPUT  VARIABLE  P-ARRAY,  ERROR  QUANTITIES 
MEASURING  RESPECTIVELY  THE  STATE  OF  CONVERGENCE  OF 
THE  P  SIMULTANEOUS  ITERATION  VECTORS.  EACH  ELEMENT  OF 
THE  ARRAY  F  IS  INITIALLY  SET  BY  SIMITZ  TO  THE  VALUE  4.0. 

EM  AS  AN  INTEGER  INPUT  VARIABLE  IS  THE  NUMBER  OF  EIGENVALUES 
TO  BE  COMPUTED,  0  .LT.  EM  .LT.  P  .LE.  N  .LE.  LDX. 

EM  AS  AN  INTEGER  OUTPUT  VARIABLE  IS  THE  NIM3ER  OF  EIGENVECTORS 
COMPUTED  THROUGH  KM  ITERATION  STEPS. 

X  AS  A  REAL  N-BY-P  INPUT  ARRAY  IS  A  SET  OF  P  OPTIONAL  INITIAL 

APPROXIMATE  EIGENVECTORS  X(I,1) .  X(I,P),  I  =  1,  .... 

N,  INTERPRETED  BY  SIMITZ  IF  KM  IS  NEGATIVE. 

X  AS  A  REAL  N-BY-P  OUTPUT  ARRAY  IS  A  SET  OF  EM  EIGENVECTORS 

X  (I  .  1 )  ,  ....  X  ( I  ,  EM)  ,  I  =  I . N ,  COMPUTED  THROUGH 

ABS (KM)  ITERATION  STEPS  WITH  THE  REMAINDER  OF  X  CONSISTING 
OF  P  -  EM  APPROXIMATE  EIGENVECTORS.  THE  P-SQUARE  MATRIX 
WHOSE  ( J ,  L)  ENTRY  IS  IP(N,  X(1,J),  X ( 1  ,L))  IS  THE  IDEN¬ 
TITY  MATRIX;  THAT  IS,  THE  EIGENVECTORS  OF  A  ARE  ORTHO¬ 
NORMAL  RELATIVE  TO  THE  INNER  PRODUCT  IP. 

LDX  IS  AN  INTEGER  INPUT  VARIABLE,  THE  EXTENT  IN  THE  LEADING 
DIMENSION  OF  THE  ARRAY  X  IN  THE  CALLING  PROGRAM. 

D  IS  A  REAL  OUTPUT  P-ARRAY  OF  WICH  D(l) . D (EM)  ARE  THE 

EIGENVALUES  OF  A  LARGEST  IN  MAGNITUDE  IM  DECREASING  ORDER 

CORRESPONDING  TO  THE  COMPUTED  EIGENVECTORS  X ( I , 1 ) . 

XU,  EM),  I  =  1 . N.  D  ( EM+ 1 )  ,  _ D(P-l)  CONTAIN 

APPROXIMATIONS  TO  PROGRESSIVELY  SMALLER  SUCH  EIGENVALUES. 

D (?)  CONTAINS  THE  MOST  RECENTLY  COMPUTED  VALUE  OF  E,  WHERE 
THE  INTERVAL  (-E,  E)  IS  THE  INTERVAL  OVER  WICH  THE 
CHEBYSHEV  ACCELERATION  WAS  PERFORMED. 

m  IS  A  REAL  INPUT  VARIABLE  UNUSED  BY  SUBROUTINE  SIMITZ.  WK  IS 
INCLUDED  ONLY  FOR  COMPATIBILITY  WITH  PREVIOUS  VERSIONS. 

WORKING  STORAGE  IS  ACCOMMODATED  BY  AUTOMATIC  ARRAYS  WITHIN 
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C  SIMITZ. 

C 

C  OTHER  PROGRAVMING  INFORMATION 
C 

C  THIS  VERSION  OF  SIMITZ  IS  DESIGNED  TO  OPERATE  IN  THE  SOFTWARE 
C  ENVIRONMENT  FURNISHED  BY  THE  SLATEC  MATHEMATICAL  SUBPROGRAM  LIBRARY, 
C  VERSION  3.1. 

0 

C  THE  PERFORMANCE  OF  SIMITZ  IS  STRONGLY  DEPENDENT  UPON  THE  CHOICE 
C  OF  INPUT  PARAMETERS  AND  UPON  THE  CAREFUL  PREPARATION  OF  THE 
C  SUBPROGRAMS  IP  AND  OP.  THE  USER  SHOULD  CONSIDER  USING  HIS  OWN 
C  ACTIVE  SUBROUTINE  INF  TO  MONITOR  PROGRESS  OF  SIMITZ  RELATIVE  TO 
C  HIS  CHOICE  OF  INPUT  PARAMETERS  IF  NO  INFORMATION  IS  OTHERWISE 
C  AVAILABLE  CONCERNING  THE  LOCATIONS  OF  THE  RELEVANT  EIGENVALUES. 

C  RECALL  THAT  SIMITZ  MAY  BE  REENTERED  WITH  KM  .LT.  0  WITHOUT  LOSS 
C  OF  INFORMATION  TO  PERMIT  CONSERVATIVE  INITIAL  CHOICES  OF 
C  ABS (KM)  ,  EPS  AND  P. 

C 

C  OTHER  PROGRAMS  REQUIRED 

C 

C  SUBROUTINE  ORTHOG  -  SUBSIDIARY  SUBPROGRAM 
C  PERFORMS  ORTHONORMALIZATION  RELATIVE  TO  THE  REAL  INNER 

C  PRODUCT  IP  BY  THE  GRAM- SCHMIDT  PROCESS  OF  THE  TERMINAL 

C  COLUMN  VECTORS  OF  THE  MATRIX  X  GIVEN  THAT  THE  INITIAL 

C  COLUMNS  ARE  IP  ORTHONORMAL. 

C  SUBROUTINE  RANDOM  -  SUBSIDIARY  SUBPROGRAM 

C  INSERTS  UNIFORMLY  DISTRIBUTED  RANDOM  REAL  VALUES  FROM  THE 

C  INTERVAL  (-1,  1)  INTO  COLUMN  VECTORS  OF  THE  MATRIX  X. 

C  SUBROUTINE  TRED2 

C  IS  THE  EISPACK  (4)  PROGRAM  WHICH  COMPUTES  A  HOUSEHOLDER 

C  TRI DIAGONAL  FORM  OF  A  REAL  SYMMETRIC  MATRIX. 

C  SUBROUTINE  IMTQL2 

C  IS  THE  EISPACK  PROGRAM  WHICH  COMPUTES  THE  EIGENVALUES  AND 

C  ORTHONORMAL  EIGENVECTORS  OF  A  SYMMETRIC  TRI DIAGONAL  MATRIX. 

C  SUBROUTINE  XERRWV 

C  PROCESSES  AN  ERROR  (DIAGNOSTIC)  MESSAGE. 

C  FUNCTION  R1MACH 

C  RETURNS  SINGLE  PRECISION  MACHINE  DEPENDENT  CONSTANTS. 

C  FUNCTION  IP 

C  IS  DESCRIBED  ABOVE. 

C  SUBROUTINE  OP 

C  IS  DESCRIBED  ABOVE. 

C  SUBROUTINE  INF 

C  IS  DESCRIBED  ABOVE. 

C 

C  METHOD 
C 

C  SIMITZ  REPRESENTS  RESULTS  OF  EXTENSIVE  MODIFICATIONS  AND  TESTS 

C  OF  SUBROUTINE  SIMITZ  (1) ,  A  FORTRAN  66  TRANSLATION  OF  THE 

C  ALGOL  60  PROCEDURE  RITZIT  (3) .  THE  BASIC  RUT I SHAUSER- RE I NSCH 

C  ALGORITHM  IS  PRESERVED. 

C  *  *  *  REFERENC  ES  (1)  PAUL  J.  NIKOLAI,  ALGORITHM  538  -  EIGENVECTORS  AND 
C  EIGENVALUES  OF  REAL  GENERALIZED  SYMMETRIC  MATRICES  BY 

C  SIMULTANEOUS  ITERATION,  ACM  TRANS.  MATH.  SOFTWARE  5 
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C  (1979).  118-125. 

C  (2)  BERESFORD  N.  PARLETT ;  THE  SYMMETRIC  EIGENVALUE 

C  PROBLEM,  PRENTICE-HALL,  ENGLEWOOD  CLIFT’S,  1980. 

C  (3)  HETNZ  RUTISHAUSER,  SIMULTANEOUS  ITERATION  METHOD  FOR 

C  SYNMETRIC  MATRICES,  NUMER.  MATH.  16(1970),  205-223. 

i 4 )  B.T.  SMITH  ET  AL ,  MATRIX  EIGENSYSTEM  ROUTINES  - 
EISPACK  GUIDE,  2-ND  ED. .LECTURE  NOTES  IN  COMPUTER 
SCIENCE  6.  SPRINGER-VERLAG,  NEW  YORK ,  1976. 
COROUTINES  CALLED  1MTQL2  .ORTHOG.Rl MACH, RANDOM, TRED2  .XERRWV 
C»**END  PROLOGUE 
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